Quantum swapping of immiscible Bose-Einstein condensates as an alternative to the 

Rayleigh- Taylor instability 
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We consider a two-component Bose-Einstein condensate in a quasi-one-dimensional liarmonic trap, wliere 
the immiscible components are pressed against each other by an external magnetic force. The zero-temperature 
non-stationary Gross-Pitaevskii equations are solved numerically; analytical models are developed for the key 
steps in the process. We demonstrate that if the magnetic force is strong enough, then the condensates may swap 
their places in the trap due to dynamic quantum interpenetration of the nonlinear matter waves. The swapping 
is accompanied by development of a modulational instability leading to quasi-turbulent excitations. Unlike the 
multidimensional Rayleigh-Taylor instability in a similar geometry of two-component quantum fluid systems, 
quantum interpenetration has no classical analogue. A crossover between the Rayleigh-Taylor instability and 
the quantum interpenetration in a two-dimensional geometry is demonstrated. 
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I. INTRODUCTION 

Bose-Einstein condensates (BECs) of dilute ultracold gases 
are ideal systems to study a wide range of many-body quan- 
tum phenomena both experimentally and theoretically. Re- 
cently, there has been much interest in the quantum hydro- 
dynamics of two-component BECs 1 1-5] and, in particular, in 
the development of quasi-hydrodynamic instabilities at the in- 
terface separating the immiscible components |0-[T^ . Among 
these phenomena, the quantum Rayleigh-Taylor instability 
(RTI) is, probably, one of the most representative and fasci- 
nating (HKl^l- The RTI releases the excessive potential en- 
ergy stored in the system and transforms it into kinetic en- 
ergy of the flow, by means of interface waves that grow into 
mushroom-shaped bubbles. In the classical case, the RTI de- 
velops when a light fluid (gas, plasma) supports a heavy one 
in a gravitational field, real or effective, as is the case, for ex- 
ample, in inertial confined fusion and Supernovae lfl3l - [T6ll . In 
classical hydrodynamics, two immiscible fluids with negligi- 
ble diffusion (e.g., water and oil) may exchange places and 
reduce the potential energy only with the help of the RTI. 

The quantum counterpart of the RTI was suggested by 
Sasaki et al. Jst] for a system of two phase-separated BECs 
with opposite projections of the hyperfine spin placed in a 
non-uniform external magnetic field, which pushes the com- 
ponents against each other. Recent results on the linear and 
nonlinear stages of the quantum RTI may be found in Refs. 
lUlll^l- Still, contrary to the classical case, the RTI is not the 
only way for the immiscible BEC components to exchange 
places. As we show in the present paper, the transformation of 
potential energy into kinetic may also happen by way of quan- 
tum interpenetration of the immiscible condensates. A similar 
effect was encountered recently in Ref. [17] in a simulation 
of a system consisting of a layer of ^^Rb sandwiched between 
two layers of '^^Rb in a 2D pancake-shaped trap. Upon abrupt 
increase of the scattering length of the **^Rb component, it 
expanded from the center to the outer parts of the trap. The 
resulting evolution of the density pattern did not resemble the 
RTI at all. Although the process was presented as develop- 
ment of the RTI, we would like to point out that Ref. iflTIl 
actually demonstrated dynamical quantum interpenetration of 



the BEC components moving from the unstable to stable po- 
sitions in an almost one-dimensional (ID) regime. 

The purpose of the present paper is to clarify the purely 
quantum nature of the encountered phenomenon. We demon- 
strate the possibility of quantum dynamical interpenetration 
of two immiscible BECs pressed against each other by means 
of an external potential in ID geometry, which eliminates 
the intrinsically multidimensional RTI. We also demonstrate a 
crossover between the quantum dynamic interpenetration and 
the RTI in the two-component BEC system for a 2D case. 

The paper is organized as follows. In Sec. |II]we formu- 
late the basic equations and parameters describing the prob- 
lem. In Sec. [nil we present numerical results for the ID case, 
accompanied by simplified analytical models to advance the 
understanding. In the process of doing this, we find an expres- 
sion for the breathing mode frequency of two phase-separated 
BECs, a result that has so far been missing in the literature. In 
Sec. |IV] we reproduce the effect in a 2D geometry and com- 
pare the process of quantum interpenetration to the quantum 
RTI. Finally, in Sec.|V] we summarize and conclude. 



II. BASIC EQUATIONS DESCRIBING A 
TWO-COMPONENT BEC 



We consider a harmonically trapped ferromagnetic spin-1 
BEC with equal population of the |1,±1) components and 
suppressed population of the |1,0) component, where \F,mf) 
is the hyperfine state with full magnetic moment of the atom 
F and its projection nif. The possibility of experimental re- 
alization of such a system has been demonstrated, e.g., in 
Ref. 12] . Since the spin exchange process — 1)) 

(|1,0),|1,0)) is suppressed, the BEC becomes effectively 
two-component. For the most part of the paper, we will be 
interested in the ID case when the BEC is tightly confined in 
the X and y directions by a harmonic potential; a 2D geometry 
will also be studied for comparison. The two components dif- 
fer only by the quantum number nif, so that equal population 
Ni = N2 = N of the components results in a symmetric ground 
state of the system with respect to the middle-plane z = 0. 

The BEC is subject to the linear Stern-Gerlach poten- 



2 



tial produced by an external magnetic field gradient, which 
presses the components against each other. In the mean-field 
approximation the wave functions of the binary BEC obey the 
Gross-Pitaevskii (GP) equations: 
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where the tildes indicate dimensional variables. The atom 
mass m is the same for both components. The total exter- 
nal potential V/(f,z) provides confinement in the trap along 
the X, y, and z axes with frequencies O)^, O),,, and (£),, respec- 
tively, and takes into account the magnetic Stern-Gerlach po- 
tential pressing the BEC components against each other In 
the following, for the most part only the frequency (£>, will be 
involved explicitly in our calculations. We assume tight con- 
finement along the y axis, o,, ^ O,, which produces a Gaus- 
sian density profile in the y direction ifTsll . In the x direction, 
we consider either tight confinement, cOx S> (O,, which leads 
to ID dynamics of the system, or — 0, which is needed to 
study a 2D geometry. (The exception is the 2D calculation 
presented in Fig. 15, where a finite (Ox was used.) As a result, 
the potential term Vj in Eqs. ([U, ^ reads 



Vj{t,z) = m(olf/2 + {-iy^BB'{t)z/2, 



(3) 



where Hb is the Bohr magneton and j = 1,2 is the number 
used to label the BEC component. The Stern-Gerlach poten- 
tial is turned on at the initial time instant f = 0, B'(f ) = B'Q{t), 
where B' = const indicates the gradient of the magnetic field 
magnitude, and 9 (I) is the Heaviside step function. In the ID 
case, the interaction parameters gij are related to the respec- 
tive i-wave scattering lengths atj as gij = Ant'?' aij / {naxaym) , 
where Qx and Uy are the oscillator length scales in the x 
and y (tight) directions of the trap ifisll . = h/mcOx, a^. — 
h/m(Oy. In the 2D geometry the interaction parameters are 
gij = AnfP'aij/ {nuym). 

We assume that the inequality g\2 > gngii is satisfied, 
which is the condition for phase separation. With g\i = g22 
and A^i ~ Ni, the hydrostatic equilibrium for B' ~ Q implies 
that the densities of both BEC components are equal in the 
bulks, and specifically, noi — no2 = no, where no is the peak 
bulk number density close to the interface in the Thomas- 
Fermi approximation (TEA) with neglected quantum pressure 
(quantum surface tension). Initially the components 1 and 2 
are located at z < and z > 0, respectively, so that a positive 
B' in Eqs. ([T]i, (|2|i pushes the condensates against each other. 

We introduce dimensionless variables as z — z/a^, t = 
f a),/2, and i/y = y/^/ ,yno; then Eqs. ([B, ^ read 
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dU, (|5]l contain only three characteristic parameters of the sys- 
tem: 
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The parameter 7 measures the relative repulsion between 
BECs, b is the dimensionless external magnetic force acting 
on the BECs at f > 0, and Rq is the TEA size of each con- 
densate in the ground state at f < 0. The TEA for each BEC 
component requires two conditions: (i) The system size in the 
z direction is much larger than the respective oscillator length 
scale a,; (ii) the interface width is much smaller than the size 
of the system. The parameter Rq also determines the healing 
length in the center of the trap fl8*l, |o = h/y/2mgno = u^/Rq, 
which is expressed in dimensionless units as 



R. 



(9) 



The characteristic penetration depth of BEC density profiles 
may be interpreted as the interface width; it is estimated as 
^int = |o/v^ lll3>[3l, or, in dimensionless units 



Thus, the TEA holds under conditions 

/?o»l, |^-V7^g»l. 

'^int 



(10) 



(11) 



The dimensionless speed of sound c, at the peak density no is 
given as — \/2Rq. 

At f < the external force b is zero, and the system is in the 
ground state which is, in the TEA, 



«1(TF)W = |V^l(z,f <0)P = 1- 
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where Vj{t,z) = z} + (— l)'fe(r). The dimensionless GP Eqs. 



for condensate 1 in the domain —Rq < z < 0, and zero oth- 
erwise; and we find a symmetric density profile for conden- 
sate 2 located in the region < z < Rq- Close to the inter- 
face, where the TEA breaks down, the interface profile can be 
approximated with the expression for a non-trapped system 

ill [US, 

"^■^"«^^'<''^^^l+exp[(-lV+i2V7/^oz]- ^''^ 



m. POSITION SWAPPING IN A ID TRAP BY DYNAMIC 
QUANTUM INTERPENETRATION 

We solve Eqs. Q numerically for different values of 
b and Rq. In our calculations we use 7= 0.01, which corre- 
sponds to the case of ^^Rb atoms in spin states |1,±1). As 
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an example of the numerical solution, Fig. 1 presents snap- 
shots of the density evolution, ni,2(f,z), for b=l and Rq = 20 
at time instants between f = and f = 450. Using the stan- 
dard propagation technique of the GP equations in imaginary 
time, we find the ground state for f < with the components 
separated in space as shown in Fig. 1(a). The TFA condi- 
tions, Eq. ( fTTT i. are satisfied for Rq ~ 20, and the numerical 
solution is in agreement with Eq. ( fT2b . except for the inter- 
face region close to z = and the BEC edges at z = ±Ro, 
where quantum pressure plays an important role. At f = 
we abruptly turn on the driving external potential [the second 
term on the right-hand side of Eq. (O], which presses the BEC 
components against each other For f > we solve Eqs. Q, 
([5j by the fourth-order Runge-Kutta method; snapshots of the 
evolution are shown in Fig. 1 (b)-(e). We also plot the time- 
resolved dynamics of the density ni (f ,z) for component 1 in 
Fig. 2; the second component is situated symmetrically, i.e. 
niitiZ) —ni{t, — z) (note, however, that at long times, numer- 
ical noise is amplified which breaks the detailed symmetry). 
The sudden turn-on of the driving field produces bulk oscilla- 
tions of the BEC components, which may be seen in Fig. 2 (a) 
at < f < 10 at the free edge of the condensate. During the 
initial oscillation stage the interpenetration is minor as shown 
in Fig. 1 (b) and the shape of the BEC components is simi- 
lar to that shown in Fig. 1 (a) with the peak density and the 
condensate width oscillating in time. 

The process of active quantum swapping starts with de- 
velopment of the first soliton of condensate 1 in the bulk of 
condensate 2 and vice versa at f w 10, as presented in Fig. 
1 (c). Eventually, more and more solitons of one BEC compo- 
nent penetrate the other, thus producing a random, seemingly 
turbulent, wave pattern. A characteristic look of this quasi- 
turbulent wave pattern is shown in Fig. 1 (d) for the time 
instant t = 40, when the centers of mass of both components 
meet at the middle plane z = 0. Interpenetration of the compo- 
nents in the form of numerous solitons may be also observed 
in Fig. 2. The driving force pushes the BEC components fur- 
ther, and they tend to separate again in the new stable positions 
z > for BEC 1 and z < for BEC 2 as we can see for a rela- 
tively long time t = 450 shown in Fig. 1 (e), and in Fig. 2 (b). 
Still, complete separation of the condensates is not observed 
even for long time intervals since the initial excess potential 
energy cannot disappear without energy loss processes, and it 
remains stored in the system in the form of kinetic energy and 
quantum pressure. For comparison. Fig. 3 presents the ground 
state of the system for the same conditions as in Fig. 1 (b)-(e), 
i.e., with both condensates placed in the stable positions of 
minimal energy from the very beginning. 

In this paper we confine ourselves to the mean-field GP 
model, avoiding thermalization effects, which means that the 
soliton size in Fig. 1 (d)-(e) has to exceed the healing length. 
For this reason, we consider only moderate strengths of the 
magnetic field b. In particular. Fig. 4 presents the center- 
of-mass coordinates of the BEC components versus time for 
different values of the magnetic field b = 0.8 — 1.5 for the 
system size Rq — 20. Surprisingly, even these moderate vari- 
ations of the magnetic field lead to dramatic changes in the 
process of quantum interpenetration. In the cases b = 1 and 
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FIG. 1 : (Color online) Snapshots of density of the BEC components 
at f = 0, 6,9.58,40,450 in dimensionless units. The inset in (a) shows 
a magnified view of the interface. Insert in (b) presents a magnified 
view of the interference pattern close to the trap edge. 
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FIG. 3; (Color online) Ground state of the system for Rq = 20 in the 
external potential b = I. 
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FIG. 2: (Color online) Evolution in time of the density ni of compo- 
nent 1 fori^o = 20 in driving field b=\. (a) Short-time development; 
(b) long-time development. Component 2 is situated symmetrically 
with respect to z = 0; however, at long time the detailed symmetry is 
spontaneously broken. 



b =1.5 the interpenetration starts rather fast and develops ac- 
tively almost from the very beginning. Taking a slightly lower 
magnetic field, b = 0.9, we observe a long preliminary stage in 
the system development, for which only bulk oscillations are 
observed with negligible interpenetration. It requires a rather 
long time interval in that case, t w 80, before active quantum 
swapping starts. Taking an even lower magnetic field, b = 0.8, 
we do not observe swapping at all, which suggests the possi- 
bility of a critical magnetic field needed to drive the process. 

After this overview, we turn in the following subsections to 
discuss the different stages in the process of quantum inter- 
penetration in detail. 



A. Preliminary stage: Bulk oscillations 

In this subsection we develop an analytical model for the 
ID bulk oscillations of the system of two immiscible BECs 
and calculate the oscillation frequency. Here we assume that 
the driving force b is sufficiently weak and does not lead to 



FIG. 4: (Color online) Positions Z\ (t) and Z2(f) of center-of-mass of 
the BEC components for Rq = 20 and * = 0.8, 0.9, 1,1.5. 



dynamic interpenetration; as we have seen, this is the case 
at least for short times if b is small enough. We first derive 
equations of motion for the center of mass (cm.) coordinates. 
Using the density-phase representation of the wave-function 
ij/j = y^exp(/^j), we reduce the GP Eqs. ([B, © to the 
hydrodynamic equations for the quantum fluid 



dffij + d~{hjVj) = 0, 
1 
2 



dfVj + ^d;v] = -^dmj - (1 + r) -dm-j 



2rrP- 



1 ~ ^" / , 1 -.9 r- 



(14) 



(15) 



where Vj = hrrC^di^j is the velocity of the quantum fluid. 
The boundary conditions for Eqs. (fT4] i. ( fTSl l demand that the 
condensates are localized. 



n^(f,±oo)=0, 
d^nj{t,±oo) = 0. 



(16) 
(17) 



The cm. coordinate of the j-th component is defined as 

Zj {t) =Nf j dzzfij (f , z) , (18) 
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where Nj = f hjdz = noa,J rijdz = const is the particle num- 
ber Using Eq.(fT4li and integrating by parts, we obtain 



dt 



dzHjVj. 



(19) 



We take the second time derivative and find 



dz{n\d^y\+ gnn\dm2) , (20) 
dz (iiidzVi + gnnidziii) . (21) 



In the above equations, the quantum pressure term provides 
zero contribution by virtue of the boundary condition, Eqs. 
(fT6] l. ( [TtT i. With the external potential Vj{t,z) given by Eq.©, 
and taking into account conservation of the particle number 
Nj, we arrive at the dimensionless form of Eqs. ( l20l i. (I2TI 1 as 



-a,?z,+z, = (-i)^-+' 



(22) 



where Fi2(f) is the dimensionless internal drag force between 
the condensates. 

Thus, the cm. dynamics of each component is governed by 
the external potentials [the first term on the right-hand-side of 
Eq. (|22]|1. and by the internal resistance (drag) force of inter- 
action between the components Fi2{t), which reads 



Fi2(f) = -(l + 7)/?o I dzn2d,m. 



(23) 



We stress that the derived equations are quite general, and 
they describe the system dynamics for all stages from bulk 
oscillations to active quantum interpenetration. In addition, 
the generalization to higher dimensions is straightforward. In 
particular, the equations are applicable to the steady state with 
no external magnetic field applied to the system, as we now 
show. 

In the TEA in ID geometry, the particle number and the 
initial cm. coordinates are given by 



Nj = iiQa 



2Ro 



%'(TF) 



:Zj{t <0) = (-l)V 



3Ro 



(24) 
(25) 



Eq. (l22l l provides a convenient way to go beyond this approxi- 
mation. In order to find the steady-state solution Zq/ = Zj{t < 
0), one may calculate the interaction force Fi2{t < 0) using 
the equilibrium profiles Eq. ( fTsT l. However, it is enough to 
note that to first order in 7, and neglecting the curvature of 
the trapping potential, the components add up to form a flat 
interface: n2(mt)(^) — "0 ^ ni(int)(z)- This immediately gives 
Fnit) = —{'i /^){^ +7)^ono' ^i^^ a result we obtain for the 
cm. coordinates 



Zo, = (-l)V.(l + 7) 



3Rq 



(26) 



In this way, we have found the first-order correction beyond 
the TEA Eq. (|25]l. 



Next, we consider an ID solution to Eqs. ( l20b . ( |2TI ) in the 
form of small perturbations of the equilibrium state. The per- 
turbations are produced by the magnetic force turned on at 
t = 0, assuming that the force is sufficiently weak, and the 
system oscillates freely at f > 0. The density profile of each 
component is represented by a bulk and an interface part. 



«l(TF)fc07 -Roit) <Z< -Aint, 

0, otherwise, 



(27) 



and correspondingly for condensate 2. Let us denote the peak 
density close to the trap center by fic{t), so that n^- — no before 
the magnetic pulse for f < 0. Within the TEA, we describe the 
unperturbed bulk profile for f < as 



"iSF)="i(TF)(f <0,Z) 



no 



1 



f2 



(28) 



We now introduce a small deviation parameter Tj (f ) for the 
central density, 



no 



(29) 



Taking into account A^i = N2= N, and assuming that the spa- 
tial dependence of n^^xF) remains parabolic, we obtain the 
bulk density profile 



"/(TF) (f >0,Z) 



(30) 



Here we introduced the time-dependent TEA radius of the 
cloud, RTFit) = 3N/[2nc{t)]. Using Eq. ^ we obtain 
Zj{t) = (-l)-''9iV/[16nc(f)], and to first order in T], 



Z,(f>0)-Z,o = -T7(f)Z 



(31) 



As above, the interface profile at f < is taken to be the 
infinite-system solution n^^jj^jj(z), Eq. ( fT3l l. and the interface 
profile at f > can be found by expanding Eq. (fT3T l to first 
order in the deviation tj of the central density. Again, we wiU 
be able to exploit the symmetry of the interface profiles, as 
we shall see shortly. Now, neglecting the term linear in the 
external field b, Eq. ( l20l i is put on the form 



niN (5| 



8i2jdzinidm2), (f>0). (32) 



Now insert the Ansatz, Eq. (|3TT l into the left-hand side of Eq. 
( |32] i. Eor the right-hand side, we have Fnit) = -(3/8)(l + 
7)7?()nQ(l + 77)^, and to first order in 77 we obtain 



[a,|+(3+27)«2]77=0. (f>0) 



(33) 



Thus, the cm. coordinates, as well as the TEA radius and the 
central density, oscillate with frequency ^/T'- 
note that the limit 7 



-27©;. Let us 
in Eq. (|33] ) cannot be taken, because 



when 7 < ^, one of the assumption made for derivation of 
Eq. ( I33] ), namely the second TEA condition in Eq. (fTTT i. fails. 
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FIG. 6: (Color online) Snapshot of the densities at time t = 2 for 
FIG. 5: The cm. coordinate Zi {t) for R^, = 20 and different values ^ and b = 5. The inset magnifies the interference pattern, 

of b, showing the bulk oscillation. 



In order to describe a smooth transition between the Hmiting 
case 7 = which is equivalent to the 1 -component BEC in a 
harmonic trap where the lowest monopole frequency is \/2ci),, 
and our case restricted by Eq. ( fTTl i. one should use proper 
solutions for the density profiles of the ground state instead of 
Eq. (O. 

In the dimensionless units the frequency is 2^/3.Q2 « 3.48, 
which may be compared with the numerical simulations: In 
Fig. 5 we present the function Zi(f) found numerically for 
Rq = 20 and different values of b, which is a magnified view 
of the region in Fig. 4 marked by the dashed rectangle. From 
Fig. 5 we find that the oscillation frequency is ss 3.5, in good 
agreement with the analytical result. According to the nu- 
merical simulations, the frequency of bulk oscillations is al- 
most independent of the magnetic field strength b. To finish 
this subsection, we point out that the bulk oscillations corre- 
spond to the counter-phase mode of the BEC interface pertur- 
bations identified in the linear theory of Ref. lfT2ll . but these 
two modes are not identical. We remind that the counter-phase 
mode involves oscillations of the interpenetration depth of the 
BEC components locally at the interface. For comparison, the 
in-phase mode corresponds to the RTI and implies bending of 
the interface as a whole with minor modifications of the inter- 
nal structure. 



B. Initial stage: onset of dynamic quantum swapping 

The main purpose of the present work is to study quantum 
swapping of two immiscible BECs pressed against each other. 
Figure 1 demonstrates the key elements of the process, and 
Fig. 4 suggests that quantum swapping may occur only for a 
sufficiently strong driving force. Still, a certain minor inter- 
penetration of the BEC components takes place for any mag- 
nitude of the driving force, even a low one, during the bulk 
oscillations. Visually, this minor initial interpenetration looks 
similar to that shown in Fig. 1 (b). As a result of this minor 
interpenetration, a small-amplitude matter wave of one com- 
ponent is injected into the "foreign" bulk of the other com- 



ponent, gets reflected from the trap edge, and forms an inter- 
ference pattern. A magnified view of the interference pattern 
close to the trap edge is shown in the inset of Fig. 1 (b) for 
b = I, R{) = 20. The pattern becomes even more pronounced 
for a stronger driving force; as an illustration, we present a 
snapshot of the system evolution for b^5 and Rq = 20 in Fig. 
6. In the process of bulk oscillations, the amplitude of the 
interfering wave oscillates as well, but these oscillations do 
not necessarily imply the onset of quantum swapping, since 
the swapping process requires unstable growth of the matter- 
wave amplitude. 

We develop now a simplified model in order to capture 
qualitatively the onset of the swapping process. The model 
has two main ingredients: A bulk part approximated by a 
TEA profile, and a tail extending into the foreign component, 
whose growth we want to study: 



ni(TF)fcf), --Ro(f)<^<0, 
ni{z,t)={ nni^){z,t), 0<z<Ro{t), 



(34) 



0, 



otherwise. 



Condensate 2 is given by a mirror image of the above, «2 (z) = 
ni(— z). In the spirit of keeping things simple, we assume the 
tail to be a constant function extending across the range of the 
foreign component: ni(inj)(z,0 = noC(0- More elaborate ex- 
pressions will not alter the qualitative results. (Note, in partic- 
ular, that a spatially oscillating part should be superimposed 
onto the ansatz, but such oscillations will cancel out in the 
subsequent spatial integrations.) 

The algebra carried out in Eqs. (fT4ll-(l22]l can be repeated 
with the difference that the integration is done only over the 
domain Q < z < Ro', the resulting equation reads 



md?f[ t/zzni(inj) 



Ra 



dzni(inj)di{Vi +gl2n2(TF)) 



(35) 



Note that the trap potential cancels part of the last term on the 
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RHS by virtue of the TFA, 



1 



(36) 



so that the equation of motion reads 
liBB'i 



Ro 







Next, note that particle conservation implies 



ri2 



'n2(TF) 



«c(f) ^ 



2/3 



no- 



(37) 



(38) 



Performing the integrations and using Eq. ( 1381 ). we end up 
with the equation of motion for the amplitude of injected mat- 
ter, 



4?C = 



4^ .A 3^ 



2/i 



(39) 



This equation of motion predicts that, when b/Ro^ y, the in- 
jected population increases exponentially on a time scale pro- 
portional to (Ib/Ro)^^^. Moreover, it predicts a critical value 
of the driving force for the onset of swapping. 
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FIG. 7: Length scale of the standing wave at the trap edge, as func- 
tion of magnitude of the driving force b. The curves present the 
analytical formula Eq. i42\ for Ro = 1 5 , 20, 30; the markers show the 
numerical result. 



Figure 7 shows good agreement with the numerical data for 
the density pattern close to the trap edge at the initial stage of 
the quantum interpenetration for = 15,20,30 and different 
values of the magnetic field. 



Developed stage: formation of solitons and quasi-turbulent 
mixing of the components 



(40) 



Close to the critical field, the balance of the two terms 
within square brackets is decisive, and the increase is super- 
exponential, in qualitative agreement with the cases b = 
0.9, 1 .0 in Fig. 4. Whether this critical field exists or whether it 
is an artifact of the crude approximation cannot be unambigu- 
ously inferred from the numerical results, but the behavior of 
the cm. coordinate in Fig. 4 is suggestive. In particular, the 
developed model predicts the critical force Zjci- « 0.1 for the 
parameters Ro = 20, 7 = 0.01 used in the simulations of Fig. 
1-3. This value is quite different from b « 0.8, apparent from 
Fig. 4. Still, we stress that the developed model is only qual- 
itative and it may not be used for quantitative predictions of 
the system dynamics. In particular, the role of bulk oscilla- 
tions were neglected in this model. 

At the end of this subsection, we obtain the characteristic 
length scale of the interference pattern. For that purpose we 
evaluate the dimensionless velocity q of the expanding matter 
wave at the trap edge from the energy conservation as 



1' 



\/2bR^, 



(41) 



where 2bRo is the excess potential energy released in the ex- 
pansion process. Then, the spatial dependence of the expand- 
ing matter wave close to the trap edge may be presented as 
exp{ztiqz), thus producing the standing wave oc exp{±i2qz) 
with wave number 2q and wavelength Asw = T^/q, 



V2bR^' 



(42) 



The process of active quantum swapping starts, when the 
first ID "droplet" appears, i.e. bright soliton of a large am- 
plitude ni ~ 1 of the component 1 develops in the bulk of the 
foreign component 2 (and vice versa). Due to the pressure 
balance, any bright soliton of one BEC component implies a 
dark soliton of the other component positioned in the same 
place. For driving field strength b — I the first soliton devel- 
ops at f « 10, as one can see in Figs. 1 (c), 2. In the following, 
the soliton evolves in an essentially nonlinear way in the bulk 
of the foreign component, interacting with the "background" 
matter, and transforming into new solitons. Still, it is possi- 
ble to trace positions of individual solitons versus time during 
rather long time intervals in Fig. 2 (i.e. the sinusoidal "tra- 
jectories" in z — f coordinates). When the number of solitons 
becomes sufficiently large, the system enters the developed 
stage of the dynamic interpenetration process. 

The developed stage of quantum swapping exhibits com- 
plicated dynamics with large velocity fluctuations. A snap- 
shot of the system evolution at f = 40, presented in Fig. 1 
(d), shows an example of density distribution at that stage. At 
this stage, we argue that the characteristic scale of the den- 
sity fluctuations must be the same as that given by the MI for 
two counter-flowing condensates. The MI was analyzed in 
Ref. 1 1 1] in the context of two partially overlapping, counter- 
flowing condensates (cf. |21]). Two initially homogeneous 
and coexisting condensates with relative wavevector q were 
found to obey the dispersion relation 
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FIG. 8: Dispersion relation for the MI for /?o = 20, fc= 1, 1.5,2 (from 
bottom to top). Solid lines denote imaginary part of the frequency; 
dashed lines denote real part. Units are dimensionless. 
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: kc/V2, which happens to be exact atq = 0; we obtain 
2n 



'^Ml - 
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(45) 



We see that at <7 = 0, Ami is equal within a constant factor to 
the interface width Ajnt = (y^/^o) ' ■ The characteristic wave- 
length, with Eq. (|4TJ inserted for q, is plotted in Fig. 9 to- 
gether with the numerically computed wavelength, calculated 
as the inverse of the number of oscillations per unit length; the 
error bars indicate the fluctuation in this quantity over time. 
We note that the MI length scale is somewhat larger than the 
length scale of the interference pattern at the initial stage of 
the quantum swapping process shown in Fig. 7, albeit compa- 
rable by order of magnitude. 

Eventually, the components tend to separate in space again 
and to occupy the positions with minimal potential energy in 
the trap as shown in Fig. 1 (e) for f = 450. In Fig. 1 (e) the 
main part of component 1 has moved from the domain of z < 
to z > 0. This tendency may also be observed in Fig. 2 (b) for 
the density evolution of component 1, rii {t,z), presented in the 
space-time coordinates on large time scales. Alt = 450 shown 
in Fig. 1 (e), separation of the condensates is not complete; we 
can see remnant excitations in the condensates created by the 
release of the excess potential energy and by the MI, which 
means that a considerable fraction of the potential energy is 
transformed into kinetic energy. In the next subsection we 
will look at energy scales in more depth. 

D. Rate of quantum swapping and energy balance in the 
process 



FIG. 9: Length scale of the multi-domain structure formed during 
the swapping process, as function of magnitude of the driving force 
b. The curves correspond to the analytical formula Eq. l l45t for Rq = 
20, 30; the markers show the average numerical values obtained from 
the simulations for Rq — 20; the error bars indicate the uncertainty in 
these values as the waves develop over a short duration of time. Units 
are dimensionless. 



The dispersion relation Eq. ( l43T l is plotted numerically in Fig. 
8 for the system size Rq = 20 and the magnetic field strength 
b = 1, 1.5,2. For all parameters, we observe an unstable do- 
main for long wavelengths, k < kc, and a stable domain at 
sufficiently short wavelengths, k > kc, where kc is the MI cut- 
off. 
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(44) 

The characteristic length scale of the instability pattern is 
determined by the wavelength Ami — 27r/fcmax correspond- 
ing to the maximal growth rate. Though it is difficult to 
find a concise analytical formula for ^max. it may be evalu- 
ated with reasonable accuracy from the numerical solution as 



The overall rate of the swapping process may be character- 
ized by the time scale. To, required for the center-of-mass of 
the BEC components to reach the middle of the trap, z — 0. 
The characteristic time Tq is plotted in Fig. 10 as a function of 
the magnetic field strength b for different sizes of the system 
^5 < Rq < 20. It is seen that reducing the magnetic field b 
or increasing the system size Rq will increase the time To re- 
quired for the swapping process. For example, by increasing 
the magnetic field strength from 1 .0 to 2.5 for Rq = 20, we de- 
crease the swapping time by almost two orders of magnitude. 
Included in Fig. 10 are data points for those parameter sets 
where the swapping process is dominated by quasi-turbulent 
mixing. As discussed in connection with Fig. 4 and Eqs. ( [39] l- 
(gOl), when Rq = 20 and b = 0.8, no swapping was observed 
for the duration of the simulation, supporting our point that 
there exists a critical field below which there is no swapping. 
The case b = 0.9, close to the critical field, appears to be gov- 
erned by two quite different sorts of process; an initial slowed- 
down growth of injected matter in addition to the subsequent 
quasi-turbulent mixing. 

To complete the study of the ID quantum swapping, we 
discuss the balance and transformation of energies during the 
process. The total energy Etot can be written as a sum of six 
contributions. 



■ + Epi2 +Et+Ei + E]^ + E, 



qp- 



(46) 



9 




FIG. 10: Time Xq after which the cm. coordinates of the two compo- 
nents meet at the midpoint, as a function of the external force b for 
different system sizes Rq = 15 — 20. 



with the pressure energy of intra-species interactions 



.7=1,2 ^ 

the pressure energy of inter-species interactions 

Epi2 = J dzil+y)Rlnin2, 
the potential energy due to the trap 

j=l.2-' 

the potential energy due to the driving magnetic force 

E (-l)' [dzbznj, 
the quasi-classical kinetic energy 



(47) 



(48) 



(49) 



(50) 



= ^ / dznj{d,<j)j) , 
and the quantum pressure energy 

7=1.2 



(51) 



(52) 



Since the TFA is well satisfied, and the interaction parameter 
7 is small, |7| ^ 1, the energy terms in Eq. ( |46] | differ by or- 
ders of magnitude, and the whole process of energy transfor- 
mation splits into two groups: The first group contains large 
energies related to pressure and the trap potential, and the sec- 
ond group involves small energies describing the essence of 
the quantum swapping. 
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FIG. 11: Energy terms as functions of time, (a) Large energies. Full 
lines represent from top to bottom at the left, iitot: total energy, E^i : 
pressure energy of intra-species interactions for component 1, E-p\2: 
pressure energy of inter-species interactions, E^i : potential energy of 
component 1 due to the trap. Dashed line, almost coinciding with 
iitot, represents the sum of the large energies, (b) Small energies. 
Topmost line represents the sum of small energies; thereafter Ei\: 
potential energy of component 1 due to the driving force, E^\ : quasi- 
classical kinetic energy of component 1, E^pi : quantum pressure en- 
ergy of component 1. 



We present the energy balance for both groups separately 
in Fig. 1 1 (a) and (b), respectively. The energies from the first 
group. Fig. 1 1 (a), are about two order of magnitude larger 
than the small energies from the second group. Fig. 1 1 (b). 
The sum of the energies from the first group, E^+Ei^ £pi2, 
represented by the dashed line, can be hardly distinguished 
from the total energy, fitot, and both values are constant within 
a very good accuracy. 



The energies from the second group provide only a mi- 
nor contribution to the total energy balance, as shown in 
Fig. 1 1 (b). At the same time, the energy balance of the second 
group reflects the physics of the swapping process. As we can 
see in Fig. 1 1 (b), the excess potential energy in the magnetic 
field is released during the interpenetration, and transformed 
into quasi-classical kinetic energy and the energy of quan- 
tum pressure isqp. The released energy is divided approxi- 
mately equally between and /Sqp. A priori, it is not obvious 
that small fluctuations of the large energies of the first group 
do not influence the energy balance in the second group. How- 
ever, taking the sum of all energies of the second group, e.g. 
Edi +£ki +£^qpi> we find only minor variations of this value 
within of less than 20% for the time span shown. 
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IV. POSITION SWAPPING IN A 2D GEOMETRY: 
CROSSOVER BETWEEN THE RTI AND QUANTUM 
INTERPENETRATION 

In this section we consider a 2D system, for which the po- 
sition swapping can occur in two ways: either by the intrin- 
sically multidimensional RTI or by the quantum interpenetra- 
tion described above. The RTI is typically a powerful phe- 
nomenon at sufficiently large length scales, but it may be re- 
duced and suppressed at small scales, e^ by transport pro- 
cesses or quantum dispersion lfl4l [Tsi |23|] . In the quantum 
systems of two-component BECs, the RTI is stabilized by the 
effective surface tension for perturbations with wavelengths 
below a certain cut-off value Art- The cut-off was calculated 
in Ref. [12] using a variational Ansatz in the limit of a wide in- 
terface, and for the introduced dimensionless variables it reads 
as 




(53) 



The theoretical analysis 111211 assumes that the inner struc- 
ture of the unstable interface is not perturbed because of the 
2D bending, which holds only for a driving force of limited 
strength. On the other hand, the quantum interpenetration 
studied above happens for a relatively strong magnetic field, 
which violates the limits of the variational analysis llT2ll . As a 
result, one should expect corrections to the cut-off value, Eq. 
(l53T l. due to the interpenetration effects, which make the un- 
stable parameter domain somewhat larger. Nevertheless, we 
use Eq. ( |53] l as the characteristic length scale for which the 
RTI becomes relatively weak. 

The purpose of the present subsection is to compare the rel- 
ative roles of the RTI and quantum interpenetration at differ- 
ent length scales close to the analytical predictions for the RT 
cut-off, Art. In particular, we obtain that the RTI dominates 
at sufficiently large length scales X > V3Art, corresponding 
to the maximal RT growth rate at the linear stage lfl2ll . In that 
case the system dynamics resembles qualitatively the simula- 
tion results of Refs. [So JJJ, and, therefore, it is not presented 
here. When the system width in the transverse direction is no- 
ticeably lower than Art, then the RTI does not develop at all 
and one should naturally expect the dominating role of dy- 
namic interpenetration. The respective system evolution is 
demonstrated for a system whose extent in the x direction is 

= O.SArt in Fig. 12, which shows the density snapshots 
of the BEC component 1, ni{x,z,t), at f = 0, 3.28, 12.0, for 
Rq — 30 and b = 5. In order to simplify the analysis, here the 
condensate is confined by a ID trap in the z direction, and we 
choose periodic boundary conditions in the x direction. The 
density of the second component 2 is not presented in Fig. 12 
for brevity, but it may be reconstructed using symmetry, which 
is unbroken until the quasi-turbulent regime develops. We also 
mark that the snapshots in Fig. 12 are strongly squeezed in 
the z direction for illustrative purposes. At the early stages of 
the system evolution presented in Fig. 12 (b), we observe a 
ID pattern of stripes similar to the quantum interpenetration 
described in the previous subsection. Still, after a while, the 
ID evolution of the system breaks down because of the multi- 
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FIG. 12: (Color online) Snapshots of density of the BEC component 
1 for = 0.5Art, /?o = 30 and = 5 at f = 0, 3.28, 12.0 in a 2D 
geometry with initially flat interface between BECs. The 2 compo- 
nent is located symmetrically. Note that proportions of the axes are 
distorted. 
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FIG. 13: (Color online) Magnified view of the region marked by 
the dashed rectangle in Fig. 13 (b): The matter- wave interference 
fringes break up into droplets due to the capillary instability at f = 
3.28 (note the adjusted color shading). The two panels show density 
111 and phase of the 1 component. The 2 component is located 
symmetrically. 



dimensional effects resembling the capillary instability 124 
The capillary instability starts developing on the fringes, close 
to the system edges. In order to make the instability onset vis- 
ible. Fig. 13 presents the part of the snapshot Fig. 12 (b) 
selected by the dashed white rectangle with equal scales in 
the X and z directions. The symmetry breaking develops with 
time into a 2D quasi-turbulent pattern of irregular droplets, as 
one can see in Fig. 12 (c) for the whole trap (the snapshot is 
squeezed in the z direction) and in Fig. 14 for the central re- 
gion with equal scales in both directions. Different from the 
ID quantum interpenetration studied in the previous subsec- 
tion, the multidimensional development of the MI produces 
2D drops of one condensate penetrating the other instead of 
the ID solitons of Fig. 1 (c)-(e). The dimensionless droplet 
size in Fig. 14 is about 1, which is comparable to the soli ton 
size in Fig. 1. As we have shown, the droplet size corre- 
lates with the characteristic length scale of the MI, Eq. ( l45T l. 



11 




-40 -30 -2D -10 10 20 30 40 



FIG. 14: (Color online) Magnified view of the region marked by the 
dashed rectangle in Fig. 13 (c): a well-developed irregular pattern 
of 2D droplets at f = 12.0. The two panels show the density «i and 
phase ^1 of the 1 component. The 2 component is located symmetri- 
cally. 
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FIG. 15: (Color online) Snapshots of density ni of the BEC compo- 
nent 1 at / = 0,4, 1 1 in the dimensionless units of the paper in a 2D 
trap with transverse frequency tt)v = lOOo), for Rq = 15, b — 1 .5. The 
2 component is located symmetrically. Note that proportions of the 
axes are distorted. 



which depends on the interface width between the BEC com- 
ponents, on the characteristic velocity of the interpenetration 
"flow" and, hence, on the magnetic field strength and the sys- 
tem size. This finding is also in agreement with the theoretical 
analysis of the quantum capillary insta bility performed in Ref. 
ll24ll . According to the linear analysis ll24|] . an infinitely long 
cylindrical rod of radius R of the condensate 1 immersed in 
the condensate 2 breaks down into droplets by the instability 
with a characteristic wavelength X ~ 4.6/?. Further nonlin- 
ear evolution of the system in Ref. |24] led to formation of 
droplets with the size comparable to the initial rod diameter. 
Extrapolating these findings to the present study, we should 
expect the droplet size in our simulations to be comparable to 
the width of the initial ID stripes (solitons) and, consequently, 
to the interface width in agreement with the present numerical 
data. 

On the other hand, the capillary instability cannot break the 
transverse symmetry when the system width is smaller than 
or comparable to the width of the matter-wave solitons pro- 



FIG. 16: (Color online) Snapshots of density nj of the BEC com- 
ponent 1 for Lx = Irt, Ro = 30, i = 5, at f = 0,3.22, 12 in a 2D 
geometry with initially curved interface between BECs. The 2 com- 
ponent is located symmetrically. Note that proportions of the axes 
are distorted. 



duced by the MI. To demonstrate this we simulate the dynam- 
ics of a 2D two-component BEC trapped in both the x and 
z directions; the simulation snapshots are shown in Fig. 15 
for RQ = 15,b=1.5 and OJ^/co, = 100 at f = 0, 4, 1 1 . Figure 
15 (a) presents the initial density for condensate 1 confined in 
the left half of the cigar-shaped trap, with n2{t,z) = n i (f , — z) ; 
the snapshot (b) shows the interpenetration process at a rela- 
tively early stage, when the solitons are already formed. The 
solitons are strongly confined in the x direction in Fig. 15, 
which suppresses the capillary instability. The snapshot (c) 
corresponds to the point when the cm. positions of both com- 
ponents are close to the middle plane, z = 0. 

Finally, taking the system width Lx close to the critical 
length scale Art we obtain a crossover between the RTI and 
the quantum interpenetration shown in Figs. 16 for the whole 
system (again, the field of view is squeezed in the z direction) 
and in Fig. 17 for the central part only (with equal scales) for 
Ro = 30 and = 5 at the time instants t 0,2.84, 12.0 and 
t = 1.88, 1.92,2.84, 12.0. Though we chose the width equal 
to Art in this simulation run, the evolution at the initial stage 
reproduces qualitatively the main features of the RTI, namely, 
unstable bending of the interface, which leads to formation 
of the "mushroom" structures shown in detail in Fig. 17 (a). 
These features of the system evolution indicate that the RT 
stability limits are wider at relatively strong magnetic fields, 
e.g. for — 5, in comparison with the prediction of Eq. ( l53T l. 
We can also recognize the 2D matter-wave stripes of conden- 
sate 1 in the bulk of condensate 2 and vice versa in Fig. 16 (b) 
at 15 < |z| < 25. However, the matter-wave stripes are rather 
weak as seen on the snapshot, so that they play a minor role 
in the system development for the width = Art, while the 
RTI bending dominates at the early stages. Surprisingly, fur- 
ther evolution of the system presented in Fig. 16 (b), (c) on 
large scales, and in Fig. 17 (b)-(d) in detail, differs qualita- 
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lively from the classical RT scenario. Because of the cap- 
illary effects and quantized vortex formation, the mushroom 
cap detaches from the stem. One can see in Fig. 17 (b) that 
the vortex-antivortex pair (vortex ring in 3D) of the compo- 
nent 1, which is located on the line z > in Fig. 17 (b) enters 
the spatial domain of the cap of the same component (note 
another vortex pair on the line z < in Fig. 17 (b), which is 
just going to enter the cap). Vortex cores under the cap of 
the component 1 capture atoms of the component 2, and visa 
versa, and thus split the cap's sides into droplets, as one sees 
in Fig. 17 (b). The mushroom evolves into a collection of 
droplets. Fig. 17 (c), and eventually breaks down into a spo- 
radic pattern, and the whole system looks more and more tur- 
bulent with time [snapshots 16 (c) and 17 (d)]. Remarkably, 
the sporadic system of droplets in Fig. 17 (d) demonstrates 
no qualitative difference from that of Fig. 14, though the lat- 
ter resulted from the RTl, while the former developed from 
the quantum interpenetration. This leads us to the conclusion 
that the quantum effects of interpenetration, surface tension 
and Ml eventually dominate in the system dynamics, at least 
on the scales about Art, and produce quasi-turbulent quantum 
mixing in the system of two immiscible BECs. 
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V. SUMMARY AND CONCLUSIONS 

We have studied the one- and two-dimensional dynamics 
of a pair of immiscible Bose-Einstein condensates subject to 
a magnetic force acting to press the components against each 
other. In ID, we find that the condensates move towards their 
respective potential minima through an essentially quantum- 
mechanical interpenetration process. The swapping process 
proceeds in several stages: Initial bulk oscillations accompa- 
nied by growing tails within the foreign bulk, subsequent ap- 
pearance of solitons penetrating into the foreign bulk, and fi- 
nally a proliferation of solitons aided by a modulational insta- 
bility which results in a mixed quasi-turbulent state, in which 
the two condensates shift towards their new equilibrium po- 
sitions. The initial part of the process is seen to happen on 
a short time scale, exhibiting exponential growth, except in 
degenerate cases. Numerical and variational calculations sug- 
gest the existence of a critical value of the driving force, below 
which there is no swapping. 

In geometries intermediate between ID and 2D, we study 
the balance between the Rayleigh-Taylor instability (RTl) and 
the quantum swapping, with the former dominating at the ini- 
tial stage of the process provided the width of the system in 
the transverse direction is greater than the cut-off wavelength 
for the RTl. In the intermediate regime dominated by quan- 
tum swapping but not truly ID, the interpenetration process 
is accompanied by capillary instabilities creating a 2D quasi- 
turbulent state. 
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FIG. 17: (Color online) Snapshots of density ni of the BEC compo- 
nent 1 at f = 2.00, 2.08, 3.22, 12.0 for the same parameters as in Fig. 
16. Red circular arrows in (a) and (b) show velocity direction of the 
BEC component 1 around the quantized vortices. 
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formance Computing Center North (HPC2N). 
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